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I. INTRODUCTION 

This paper deals with the 'inverse problem' in classi- 
cal statistical mechanics. Namely we are interested in 
determining the interaction potential of a system from 
the knowledge of its radial distribution function (RDF). 
A basic result due to Henderson [1] states that if a sys- 
tem is governed by pairwise additive interactions then 
two potentials which give rise to the same RDF cannot 
be different more than a constant term. This theorem 
provides a theoretical support to the formulation of the 
inverse problem since it demonstrates the uniqueness of 
its solution. However the existence of the solution is not 
guaranteed and furthermore the theorem does not indi- 
cate a way to find it. 

Despite this general result, the solution of the inverse 
problem for a classical dense fluid turns out to be a diffi- 
cult task to achieve. This is due mainly to the fact that 
in the high density regime the RDF is hardly sensitive to 
the detailed shape of the interaction potential and is es- 
sentially determined by its repulsive part; so the inverse 
functional relationship between the RDF and the interac- 
tion potential evidences a strong dependence of the latter 
on the input RDF. In order to expect a reliable solution of 
the inverse problem not only the input RDF must be pro- 
vided with high precision but also the underlying theory 
used to formulate the inversion procedure must be very 
accurate. As stated by Reatto in [2] the accuracy of a 
satisfactory inversion scheme has be to independent both 
from the shape of the interaction potential and from the 
density of the system under inspection. If these proper- 
ties are fulfilled then the interaction potentials of differ- 
ent systems can be consistently compared, furthermore 
any dependencies of the extracted potential on the ther- 
modynamic state can be unambiguously ascribed to the 
effects of many-body interactions. 

A generally accepted scheme for the solution of the 
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inverse problem which fulfills these features is still lack- 
ing and, in the last three decades, several authors have 
proposed different approaches. A first category com- 
prises theoretically based attempts in which the inversion 
scheme is defined on the basis of an integral equation 
theory (HNC, MHNC, etc.) of the liquid state. These 
pure theoretical approaches typically rely on some ap- 
proximation and, due to the intrinsic difficulties depicted 
above, their application provides reliable results only in a 
limited set of cases. A considerable improvement in the 
accuracy of the extracted potentials has been obtained 
by recurring to simulation assisted procedures. These 
methods attempt to determine the pair interaction start- 
ing from a guessed expression of the potential which is 
iteratively modified on the basis of the discrepancy be- 
tween the simulated pair function and the experimental 
data. A first result in this direction has been proposed 
by Schommers in [3] and later on further improvements 
have been achieved by Reatto, Levesque and Weiss in [4]; 
in this paper the authors applied the predictor-corrector 
scheme, using the MHNC equation as predictor, to the 
Lennard- Jones fluid and to a model potential for alu- 
minum. The convergence of the iterative potential to the 
correct result was found and it was checked that the use 
of a less accurate predictor (for example the one proposed 
in [3]) for the definition of the trial potential could spoil 
the accuracy of the procedure. Other results belonging 
to this class of inversion procedures comprise the em- 
piric potential structure refinement (EPSR) proposed by 
Soper [5, 6] and a solution due to Lyubartsev and Laak- 
sonen [7] . The former technique performs the refinement 
of a reference potential using a perturbation term given 
by the difference between the experimental and the sim- 
ulated structure factor; the latter propose a parametric 
dependence of the potential on a set of parameters which 
are determined by solving a large system of linear equa- 
tions. 

A further approach to the inverse problem is provided 
by a family of 'stochastic' inversion methods in which 
the solution is sought as the expected value of properly 
extracted random variables (inverse Monte Carlo). In 
this simulation scheme, given the input RDF, a dynam- 
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ical evolution law is denned with the aim to build a set 
of configurations compatible with the experimental data. 
So the solution of the inverse problem is brought back 
to the determination of a suitable transition probability 
which produces a 'Monte Carlo like' dynamic. Among 
the various attempts in this direction we mention the Re- 
verse Monte Carlo (RMC) technique due to McGreevy 
and Pusztai [8] and two 'absolute minimization' meth- 
ods proposed by Cilloco in [9] and later on by da Silva, 
Svensson, Akcsson and Jonsson in [10]. Strictly speaking 
these methods do not provide a solution of the inverse 
problem since they do not allow the direct determination 
of the interaction potential, however the configurations 
produced in the inverse Monte Carlo procedure can be 
used to compute quantity of physical interest. It is worth 
mentioning that results reported in [9] represent the first 
application of the maximum entropy principle, indicating 
a possible solution based on the measurement of the three 
body correlation function. A further contribution is due 
to a technique proposed by Almarza, Lomba and Molina 
[11, 12] where a direct solution of the inverse problem 
has been obtained by performing a continuum refinement 
procedure of a trial interaction potential. 

The purpose of this paper is to present a technique for 
the solution of the inverse problem based on the max- 
imum entropy principle (ME) [13]. ME is an effective 
tool for setting up the equilibrium distribution of a sta- 
tistical system on the basis of partial knowledge and the 
corresponding estimate fulfills the remarkable property 
of being the 'maximally noncommittal with regard to 
the missing information'. So, our solution of the inverse 
problem is based on the maximization of the configura- 
tional entropy constrained by the information codified in 
the radial pair distribution function. The procedure is 
realized inside an inverse Monte Carlo scheme and the 
interaction potential emerges as the asymptotic expres- 
sion of the transition probability. 

The contents of the paper are as follows. Section II 
contains a description of our method, in section III we 
test the method in the case of a Lennard- Jones fluid and 
for a model of liquid aluminum. Finally in section IV we 
discuss our results and present some final remarks. 



II. THEORY 

A. Statistical description of a monoatomic system 

We perform a statistical analysis of a simple 
monoatomic system with the aim to define some quan- 
tities that will be of central interest later on in the pa- 
per. Particular emphasis will be given to the concepts 
of probability, likelihood, entropy and to their mutual 
relationship. 

Consider an homogenous and isotropic system com- 
posed of point-like elements with average density p. In 
the following we will refer to this system as the model. 
Given an arbitrary configuration x of the model we can 



perform a local sampling of the elements pair function 
(PF) . This means that we select a reference element and 
divide the space in spherical shells of width Sr centred 
on it up to the maximum value tm [14]; then we count 
the number of elements in each shell and we store these 
numbers in the array rij, where i — 1, J. 

We define a probability function p(x) over the config- 
uration space of the model system and collect an ensem- 
ble of s configurations extracted according to p(x). The 
global sampling of the PF over the ensemble can be com- 
puted by evaluating the rij for each configuration a and 
summing these local samplings shell by shell, that is: 

mi = XX Q) (!) 

a=l 

Assuming that the expected PF is given by a reference 
function in we can evaluate the probability associated to 
the global sampling (1). Let us focus on a fixed shell k. 
The values of obtained in two different configurations 
are uncorrelated and, admitting that the number of shells 
is large enough, the probability of finding more than one 
element in a single measurement can be neglected; so the 
shell k follows a Poisson distribution with expected value 
spk- Since there is no correlation between different shells 
we obtain: 

W = ne--^P (2) 

»=i m% - 

so the probability associated to m is given by a product 
of Poisson distributions. This formula describes an 'open 
system', which can be realized as an open subset of a 
larger one, and elements fluctuations are possible. Con- 
versely, if we are dealing with a closed system in which 
each configuration of the ensemble satisfies the further 
constraint: 

J J 

Ul = = N p (3) 

i=l i=l 

the total number of elements is conserved and the proba- 
bility (2) is reduced to a multinomial expression (see [15] 
and references therein): 

where N = sN p and qi — p.i/N p is the normalized refer- 
ence probability distribution. Equations (2) and (4) can 
be interpreted as likelihood functions C(m, p) of the ex- 
pected PF given the observed values m. If the number of 
configurations in the ensemble is very large (s> 1 which 
implies N,m>i ~> 1) we can take the logarithm of £ and 
make use of the Stirling approximation up to the linear 
order, this gives: 

ln£-p(m, /j) ~ — m i m — ~ + ^2( m i ~ s ^i) 
]n£ M (m,n) ~ ~X] miln lv ~ ( 5 ) 
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the two formula in (5) differ for a linear term which ac- 
counts for the fluctuations of elements. 

The log-likelihood equations (5) possess a nice inter- 
pretation when the number of configurations becomes in- 
finite. Let us focus on the multinomial likelihood given 
by the second line of (5); in the asymptotic limit the av- 
erage PF converges to the probability pi = m,i/N built 
over the ensemble and the likelihood can be written as: 



lim - In Cm (m, fj) 

s^oo s 



-N, 



^2 Pi ln 



Pi 



-N p D(p\\q) 



(6) 

we recognize that the log-likelihood is proportional to 
the relative entropy D(p\ \q) (Kullback-Leibler divergence 
[16]) of the ensemble distribution p with respect to the 
reference one. The relative entropy fulfills the properties 
of being positive definite and vanishing only if p = q. 
Equation (6) implies that if the global PF built over the 
model ensemble maximizes the likelihood with the ref- 
erence function p, then, asymptotically, the distribution 
p minimizes the relative entropy respect to the reference 
probability q. We will make us of this property in the 
next section. 

It is useful to rewrite equation (6) in term of radial 
distribution functions. The model RDF g(ri) and its ref- 
erence counterpart g (ri) arc defined by normalizing the 
ensemble average and the expected reference function p 
by the average value of particle per shell, respectively. So 
we have: 



1 



g{n) = lim - 



mi 



oo s 4irpr~5r 



9o{n) = 



AirprfSr 



(7) 



Plugging equation (7) in (6) and passing to the contin- 
uum limit provides an expression for the relative entropy 
that will be widely used in the following: 



K M{g\\9o)= lim ^-\n£ M (m,p) 

s-s-oo As 



(8) 



where the extra factor ^ has been inserted to avoid a 
double counting of the number of independent distances 
between pairs of elements. The same analysis can be 
repeated starting from the first line of equation (5); per- 
forming the asymptotic limit and recasting the result in 
term of RDFs gives: 



K-p(g\\g 



lim — ln Cp(m, p) 

s— too ZS 



-II 



dr 



g(r) ln 



9{r) 
9o(r) 



{g(r)-g (r)) 



(9) 



which provides the relative entropy between the RDFs 
when elements fluctuations are taken into account. 

A last comment regards the meaning of this construc- 
tion when a uniform reference distribution % = 1/J is 
employed. In this case equation (4) provides the number 
of occurrences of the global PF (1) up to a constant factor 



and the relative entropy D(p\\q) reduces to the Shannon 
entropy [17] up to an additive constant. Expressing this 
condition in terms of RDFs supplies the measurement 
of the relative entropies (8) and (9) respect to the 'non- 
informative' reference system g = 1: 



C(2) 

M 



K M (g\\i) 



C(2) 



K v (g\\l) (10) 



exploiting equations (8) and (9) we recognize that the en- 
tropies (10) exactly reproduce the two-body contribution 
to the Boltzmann entropy expansion in the canonical en- 
semble [18] and in the grand canonical ensemble [19, 20], 
respectively. 



B. Maximum entropy solution of the inverse 
problem 

We consider a monoatomic system whose interactions 
are governed by a genuine pairwise additive potential 
4>{r) and assume that for a given condition of temper- 
ature T and density p the RDF of the system g (r) is 
known. We refer to this system as the target. The inter- 
action potential of the system is supposed to be unknown, 
only the RDF is given. 

We propose a solution of the inverse problem based 
on the maximum entropy principle [13] constrained by 
the information encoded in the RDF of the target sys- 
tem. Namely we build a probability distribution p in the 
model system which fulfills the properties of maximizing 
the Shannon entropy consistently with the condition of 
vanishing relative entropy with respect to g (r). 



K(g\\g o ) = 



(11) 



where the model RDF g(r) is obtained by averaging the 
global PF (1) over an ensemble of configuration extracted 
according to p. Formally this task is achieved by com- 
puting the maximum of the functional: 



Hp} = S{p} + aK{g{p}\\g ) 
where S{p} is the Shannon entropy: 

S{p} = - ^Pn hip„ 



(12) 



(13) 



and a is a Lagrange multiplier. The stationary point 
of (12) provides the equilibrium distribution constrained 
by the target RDF and we will show that the knowledge 
of this function allows to introduce a notion of interac- 
tion potential in the model system. This quantity will 
be identified with the target potential thus providing a 
solution of the inverse problem. 



1. Low density solution 

In the low density limit the general strategy previ- 
ously described can be easily carried out. In order to 
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evaluate the stationary point of the functional (12) we 
perform an expansion of the Shannon entropy in correla- 
tion functions. Leaving aside the ideal-gas contribution 
which does not depend on the configurational degrees of 
freedom we have: 

%} = E s(n) ( 14 ) 

n>2 

Formula (14) provides an expansion of the excess entropy 
organized in powers of the density and in the low den- 
sity limit the whole series is dominated by the two-body 
contribution S l2> . 

The solution of the inverse problem is straightforward 
and proceeds in two steps. First of all we maximize the 
two-body Shannon entropy assuming that the dynamics 
in the model system is governed by an (unknown) pair- 
wise additive potential 4> m (r). For pairwisc additive in- 
teractions the configurational part of the internal energy 
can be expressed as: 

U = P - J drg(r)4> m (r) (15) 

so the ME estimate of the two-body entropy functional 
subjected to the average value of the internal energy is 
given by the stationary configuration of the functional: 

Hg} = S? + « (f / dvg{r)^ m {r) - uj (16) 

maximizing (16) and imposing the constraint (15) to- 
gether with the thermodynamic relation (3 — dS/ dU pro- 
vides the solution: 

g (r) = e-WmW (17) 

which is the ME estimate of the two-body equilibrium 
distribution for a system with pairwise interactions [21]. 
We recognize the first order contribution in the cluster 
expansion of the RDF. 

The second step is realized by imposing equation (11) 
which allows to evaluate the ME estimate of the inter- 
action potential <p m (r) constrained by the target RDF. 
Since the vanishing of the relative entropy implies the 
equality of the two RDFs we obtain: 

/3<j> m (r) = - In g (r) (18) 

which is the ME solution of the inverse problem at low 
density. 

2. High density solution: a Monte Carlo approach 

The correlators expansion of the excess entropy (14) for 
a high density system contains, apart from the two-body 
contribution, all higher order terms. Since these quanti- 
ties are unknown a direct maximization procedure of the 
excess entropy, like the one performed in the low density 



limit, is unfeasible. However, if the interaction potential 
is pairwise additive, the RDF still codifies all the infor- 
mation needed to the solution of the inverse problem. 
This is a direct consequence of the Henderson theorem 
[1]: the RDF determines the interaction potential up to a 
constant, so its knowledge sets the whole configurational 
part of the phase space distribution function and all the 
higher order terms in the entropy expansion are theoret- 
ically determined if the two-body contribution is given. 
Anyhow, since the explicit computation of these terms 
would require the knowledge of the interaction potential, 
a direct maximization procedure cannot be performed 
and a different approach has to be adopted. 

The general strategy to achieve the entropy maximiza- 
tion is to recur to a 'Monte Carlo like' (MC) approach 
in which the configuration space of the model system is 
sampled along a random path. So, as in the standard 
Metropolis-Monte Carlo (MMC) algorithm, the dynam- 
ical evolution of the system is defined by introducing a 
notion of trial configurations and a transition probability 
between neighbour states. We shall see that the stochas- 
tic nature of the MC dynamics together with a suitable 
choice of the transition probability will allow to gener- 
ate a path in the configuration space of the model sys- 
tem which maximizes the excess entropy (13) consistently 
with the relative entropy constraint (11). 

Let us define the building blocks of this procedure. As- 
sume that we have performed s MC iterations. For each 
point of the path we compute a local sampling of the PF 
and sum up these measurements in the global pair func- 
tion (1). Then we select a reference particle and compute 
a local sampling of the PF n {1) , at the same time the par- 
ticle is randomly moved and the new local sampling of 
the PF is stored in the array n (2) . This procedure pro- 
vides two different samplings of the global PF at the level 
s + 1: 

m (1) = m + n (1) m <2) = m + n (2) (19) 
the trial configuration m {2) is accepted with a probability: 

P m <D-> m C2) =min(l,/(m< 1 ',m< 2 ')) (20) 

where / is the transition probability which determines 
the stochastic evolution law. The iteration of this proce- 
dure allows to generate the whole ensemble of configura- 
tions of the model system. 

Now we impose the constraint (11). To achieve this 
task we define the transition probability by the require- 
ment that the global PF (1) built along the path maxi- 
mizes the likelihood function (5) with the reference pair 
function fi, defined in term of the target RDF via the 
relation: 

IH = Airprfg (ri)Sr (21) 

If we are able to impose this condition then equations 
(6) guarantees that, asymptotically, the relative entropy 
between the model and target RDFs vanishes and the 
constraint (11) is satisfied. For this purpose we try to 
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guess a formula for the transition probability written in 
term of a likelihood ratio: 



f = e~ 5X , where SX 



In^^iO (22 ) 



so trial samples with a likelihood higher than m (1) are au- 
tomatically accepted, otherwise they are accepted with a 
probability given by /. For s> lwe can make use of the 
Stirling approximation (5) for the log-likelihood terms in 
(22). Moreover, since the rii are of order 1 while the nii 
are of order s we can expand in power of s the logarithms 
appearing in (5). Performing this approximation to the 
first order in 1/s provides: 



to; 



(23) 



this formula can be obtained starting from both the ex- 
pressions for the log- likelihood given in (5), so the tran- 
sition probability (23) turns out to be invariant respect 
to the boundary condition imposed in the model system. 
Equation (23) computes the difference among n (1) and 



v.- 



weighting each shell with a term: 



= In 



(24) 



that represents the 'error' after s iterations between the 
reference and the measured values of the global PF. So 
8X realizes a feedback in the model system, since it be- 
haves as a controller which selects the configurations in 
the model ensemble on the basis of the error (24). This 
controller operates only by considering the error in ac- 
tual state s and, adopting the common language of the 
feedback control systems [22], we will call this quantity 
a 'proportional' controller. 

The transition probability (23), realized as a propor- 
tional controller, suffers of a difficulty which is commonly 
encountered in many feedback controlled systems when- 
ever the controller is realized only through a proportional 
term: the presence of an offset between the measured pro- 
cess variable and the target reference function. Indeed a 
MC simulation built with this transition probability pro- 
duces a model RDF which is a 'biased' reconstruction of 
the target one, so the formula guessed for 5X turns out 
to be inadeguate to enforce a complete maximization of 
the likelihood function C(m,fi). A possible solution of 
this problem can be accomplished by realizing the con- 
trol mechanism as a proportional-integral controller (PI) 
[22]. So we propose a modified expression for 5X given 
by: 



5X = J2(r 



(25) 



where Ui is a function of the error (24) which depends on 
three different contributions: a proportional term that 
determines the reaction to the current error, an integral 



term which keeps into account the sum of all the former 
ones and a background value which allows to include a 
priori knowledge on the system. The output of the PI is 
given by a weighted sum of these three quantities: 



a=l 



(26) 



where k p and kj arc the (s dependent) coefficients of the 
proportional and of the integral terms. 

A transition probability defined in term of the PI (26) 
ensures that the model RDF converges to its reference 
value. Furthermore the implementation of this controller 
allows one to define an interaction potential in the model 
system. In fact, as long as the measured PF converges 
to its reference value, the error (24) goes to zero. In 
this limit the proportional term of (26) becomes neg- 
ligible and the integral approaches to a constant finite 
value. Formally we can define the model potential as the 
asymptotic limit of PI controller (f3 = 1/fceT): 



lim u 

S— 7-OC 



(27) 



So the MC dynamics built with the PI control system 
behaves as a constructive tool for the computation of the 
model potential. During a MC simulation the model sys- 
tem is subjected to a transient dynamical phase in which 
the transition probability evolves during the path; as long 
as the path proceeds the PI builds the model potential 
(27) and the transition probability approaches to a sta- 
tionary regime. Once the equilibrium has been reached 
the system evolves according to a stationary transition 
probability and behaves as a Markov chain, in which the 
potential is given by (27). 



3. Computation of the PI coefficients 

Let us come back to the issue of the correct definition 
of the coefficients k p and fc/. Usually the PI parame- 
ters are tuned with the aim to ensure a fast and stable 
convergence of the measured process variable to its ref- 
erence value. In this case we propose a criterium, for 
fixing these parameters, which comes again from statisti- 
cal considerations. We observe that if the model system 
is sampled with the expected distribution (2), the global 
PF approaches to s/z, as long as s increases. So we in- 
troduce the reduced variables Xi defined by: 



m _ 1 | rrij - 



SfJ-i 



1 + Xi 



(28) 



and we expand the distribution function (2) in series 
around Xi = 0. Performing this expansion together with 
the usual Stirling approximation provides: 



t 



^2 



(29) 
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so for large values of s the global PF is distributed accord- 
ing to a product of J gaussian distributions [23]. Since 
the reduced variables Xi are distributed according to a 
standard normal distribution, the variable defined as: 



J , n2 

(mi - s fa) 



S^i 



(30) 



follows a %-square distribution with J degrees of freedom. 

So we define the PI coefficients in order to implement 
the condition \ 2 = 1. Enforcing this condition in the 
model system guarantees that the global PF has the cor- 
rect fluctuation around its average value and excludes 
spurious correlation among different shells. This can be 
done by introducing a new PI which performs a dynamic 
control on the coefficients k p and fcj, so we set: 



(x(.)-l)+^E(x(a)-^ 
a = l 

^^.fe-^^fe- 1 ) ( 31 ) 

a=l 

where c x , c 2 , c£i , d 2 are the PI parameters. Further details 
concerning the implementation of this control mechanism 
will be given in section III. 

III. APPLICATIONS 

In order to illustrate the features of the technique here 
proposed we have solved two systems which have been 
widely analyzed in the literature concerning the inver- 
sion methods [4, 1 1] : a simple Lennard- Jones fluid and a 
model for liquid aluminum [24]. 

We briefly describe the general strategy adopted in the 
analysis of both systems. The target RDF has been eval- 
uated recurring to a MMC simulation in the NVT en- 
semble. The configuration space of the target system is 
a cubic volume of linear length L with N p point-like par- 
ticle and the periodic boundary conditions together with 
the minimum image convention have been adopted. The 
target potential <j)(r) is truncated at L/2 and the system 
evolves starting from an FCC lattice; after about 5 x 10 2 
MMC steps the energy of the system approaches to a con- 
stant value and the system evolves around equilibrium. 
Once at equilibrium a local sampling of the PF is per- 
formed for each configuration and the average value of 
/j, is built, then the target g(r) is computed. Due to the 
minimum image convention this method provides a reli- 
able RDF up to the edge value Tm = L/2. The error on 
the target RDF can be estimated by dividing the whole 
simulation in blocks and by computing the standard de- 
viation Sg(r) between the blocks. 

Once the g(r) has been computed the inverse procedure 
for the determination of the pair potential described in 
section II can be applied. The model system is realized 
exactly as the target one, so the configuration spaces of 



the two systems are identical. The PI coefficients are 
dynamically defined by equations (31) which ensure the 
correct equilibrium fluctuation of the model RDF. A di- 
rect analysis of the system response evidences that an 
optimal choice of the parameters appearing in (31) is 
given by: 



fci s) 



k\ s) = 5 x lO- 3 /^) 



= {xU - 1) + 1 >< 10 - 3 E (4., - 1) 



a=l 



(32) 



where the ratio between k p and kj has been set to a 
constant value. This choice guarantees a smooth conver- 
gence of the measured PF to the target reference value. 
It is worth noting that performing a different choice (in- 
side a range of values which does not produce an oscillat- 
ing behavior) has only the effect of changing the rate of 
convergence of the model system but does not affect the 
convergence value. Furthermore, the same set of param- 
eters given by (32) have been used both in the analysis 
of the Lennard-Jones fluid and of the liquid aluminum, 
providing an equally good convergence independently of 
the details of the system. 

We observe that the target RDF of both the sys- 
tems under inspection exhibits a hard core structure, i.e. 
g(r) = for r < r . This information can be imposed in 
the model by introducing a hard sphere (HS) background 
potential, u (a) = oo for r < r a and zero otherwise, which 
initializes the PI controller (26). Due to this term any 
trial configuration containing particle at a distance lower 
than r is automatically rejected. Consistently with the 
background potential, we choose an equilibrium HS con- 
figuration as the starting point for the MC path. Then 
the reverse procedure starts and the system evolves ac- 
cording to the transition probability (25); after each it- 
eration we compute the output of the PI (26) and the 
expression of the transition probability is updated. Since 
the RDF of the starting configuration is noticeably dif- 
ferent from the reference value, the % 2 is sensibly higher 
than 1 and the PI coefficients (32) grow very fast; this 
phase is characterized by a highly non stationary dynam- 
ical evolution of the transition probability (25). 

In order to improve the convergence of the model po- 
tential it is convenient to split the simulation into two 
phases. So, when the \ 2 has reached a value quite close 
to 1 the actual configuration and the final expression of 
the PI output are stored in a file and we stop the simula- 
tion. Then these quantities are used as input values for 
the background potential and for the initial configuration 
and we start the 'refinement phase'. Since the system is 
closer to equilibrium, the PI (32) works in a different 
regime with respect to the previous phase; so the system 
evolves smoothly to equilibrium and the transition prob- 
ability approaches to its asymptotic value. This phase 
can be repeated many times in order to obtain a better 
refinement of the model potential. 

As a final check of the goodness of the results provided 
by this procedure we perform a standard MMC Simula- 



tion using the model potential and we compare the cor- 
responding RDF with the target one. If the difference of 
the two RDFs is not bigger than their intrinsical noise 
we conclude that the model potential (27) is equivalent 
to the target one and the reconstruction procedure stops; 
otherwise further refinement phases could be needed. 



A. Lennard-Jones potential 



The system is defined by a Lennard-Jones potential: 



<pL.j{r) = 4e 



(7)"-®' 



(33) 



with argon-like parameters a — 3.405 A and e/ksT = 
119.76. The MMC simulation for the determination of 
the target RDF is performed on a system of 864 parti- 
cle at the reduced density p* = pa 3 = 0.84 and reduced 
temperature T* = ksT/e = 0.75, near the triple point. 
The g(rt) has been evaluated up to r* = r/a = 5.0 which 
corresponds to L/2; the width of the shells for the mea- 
sure of the g{ri) was Sr = 2.5 x 10~ 2 A and the number 
of measured points was 686. We performed 2 x 10 4 cy- 
cles after equilibration. The experimental error on the 
RDF was estimated by computing the standard devia- 
tion Sg(ri) between 50 blocks of 4 x 10 2 cycles each. The 
largest value for Sg(ri) was about 2 x 10~ 2 with an aver- 
age value of 7 x 10~ 3 . 

The complete inverse simulation procedure took 2.4 x 
10 4 iteration. A first phase of 6 x 10 3 steps was performed 
starting from the FCC lattice and then the refinement 
phase was repeated three times for 6 x 10 3 steps each. The 
result for the interaction potential is reported in Fig. 1, 
the maximum difference between the model potential and 
the Lennard-Jones reference one was less than 5 x 10~ 2 
with an average value of 1 x 10 ~ 2 . The average difference 
between the model and the target RDFs was equal to 
4 x 10~ 3 ; this value is inside the average noise of the 
RDF, so the model potential of Fig. 1 can be considered 
identical to the Lennard-Jones one. 




FIG. 1. Results for the Lennard-Jones system. The target 
potential (continuous line) and the model potential (filled cir- 
cles) are plotted. 



The inverse simulation procedure took 2.6 x 10 4 itera- 
tion. A first phase of 6 x 10 3 steps was performed starting 
from the FCC lattice and then the refinement phase was 
repeated twice for 6 x 10 3 steps each and once for 8 x 10 3 
steps. The result for the interaction potential is reported 
in Fig. 2, the maximum difference between the model 
potential and the Al model reference value was less than 
2 x 10 -2 with an average value of 7 x 10~ 3 . Even in this 
case the average difference between the model and the 
target RDFs is inside the typical noise of the RDF. Ana- 
lyzing Fig. 2 we observe a difference between the target 
and model potential of the order of 1 x 10~ 2 in the range 
from 7 to 11 A. This error is due to a correlated statis- 
tical fluctuation in the reconstruction procedure and can 
be further reduced by increasing the information content 
in the target RDF used as input. 



B. Model potential of aluminum 

The system is defined by a model potential for liquid 
aluminum [24] . The MMC simulation for the determina- 
tion of the target RDF was performed on a system of 864 
particle at the density p = 0.0527 A" 3 and T = 1051 K. 
The g(ri) has been evaluated up to r — 12.70 A which 
corresponds to L/2; the width of the shells for the mea- 
sure of the g(rt) was Sr = 2.5 x 10 -2 A and the number 
of measured points was 508. We performed 2 x 10 4 cy- 
cles after equilibration. The experimental error on the 
RDF was estimated by computing the standard devia- 
tion Sg(ri) between 50 blocks of 4 x 10 2 cycles each. The 
largest value for Sg{ri) was about 2 x 10~ 2 with an aver- 
age value of 6 x 10~ 3 . 



IV. DISCUSSION AND CONCLUSIONS 

The method presented so far supplies an accurate so- 
lution to the issue of determining the interaction poten- 
tial from the radial distribution function. This technique 
bases its theoretical support on the maximum entropy 
principle of information theory which provides a general 
tool for the statistical inference on the basis of partial 
knowledge. The method is formally summarized by equa- 
tion (12) which describes the maximization of the con- 
figurational entropy (S term) constrained by the infor- 
mation codified in the target system (K term). The ME 
solution is sought inside a Monte Carlo scheme where 
the maximization of configurational entropy is realized 
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FIG. 2. Results for aluminum. The aluminum potential (con- 
tinuous line) and the model potential (filled circles) are plot- 
ted. 



through the MC random displacements and the accep- 
tance criterion for the trial configurations is built con- 
sistently with the physical input provided by the target 
RDF. The potential emerges as the asymptotic expres- 
sion of the transition probability and, for pairwise po- 
tentials, it reproduces completely the interactions of the 
target system. This method fulfills some nice properties 
that, in our opinion, make it a valid tool for the extrac- 
tion of potential. Actually the expression of the transi- 
tion probability (25) is motivated only by the constraint 
(11) and does not rely on any ulterior hypothesis con- 
cerning the physical nature of the target system, so we 
expect that the general strategy depicted in the present 
paper could be of wide applicability. Nevertheless, the 
convergence of the model potential is ensured by a feed- 
back control mechanism and the coefficients of (26) are 
tuned by an independent PI which operates a control 
on the fluctuation of the model RDF around the target 
reference value. This further controller avoids spurious 
correlations in the model RDF and guarantees that no 
information, besides the one codified in the target RDF, 
is transferred to the model during the simulation. 

Results of section III show that the extracted potential 
(27) accurately reproduces the original pair interaction 
both for the Lennard- Jones fluid and for the liquid alu- 
minum model. A comparison between these results and 
the ones presented in [4] and [11] evidences a very satis- 
factory accuracy, despite a modest computational effort. 
This level of agreement turns out to be highly remarkable 
since the systems lie in the high density region of the state 
space where it is expected that the RDF should be quite 
insensitive to the details of the interaction; moreover the 
aluminum potential exhibits well defined oscillations even 



at short distances, where f3(j)(r) is still positive. As a fur- 
ther control we have verified that the method provides 
the correct results in a different region of the (p, T) plane; 
so the procedure described in section III has been re- 
peated for a Lennard- Jones fluid at p* = 0.5, T* = 1. As 
expected, the interaction potential approaches the cor- 
rect result with a convergence rate even faster than in 
the high density case (about 1 x 10 4 steps were needed 
to obtain an accuracy comparable with the result of Fig. 
1). This analysis indicates that our procedure for the 
solution of the inverse problem provides reliable results 
independently both from the density of the system and 
the shape of the potential under inspection, so it fulfills 
the requirements of a 'satisfactory inversion scheme' as 
stated in [2]. 

The interpretation of the transition probability as a 
feedback controller represents a key point for the accom- 
plishment of the solution discussed in the present paper. 
Actually, the adoption of this point of view motivates 
the introduction of the integral term and gives rise to 
the model interaction potential (27). We want to point 
out that this is not the only way to impose the constraint 
(11). For example, the offset between the target reference 
function p and the model global PF can be made null by 
using a proportional controller with an infinite value of 
the coefficient k p . Pursuing this approach leads to the 
'pure minimization methods' [9] and [10] in which only 
trial configurations with a higher likelihood function (or 
with a lower x 2 m the language of [10]) are accepted. 
The drawback of this approach is that, due to the lack- 
ing of the integral term, the interaction potential cannot 
be directly computed. 

We conclude our discussion with some comments con- 
cerning the extension of this procedure to other systems 
than the simple monoatomic fluid analyzed in the present 
paper. The method is based on ME principle which holds 
for any system at equilibrium. For simple fluids a K term 
realized as the relative entropy (8,9) between the RDFs 
is able to constrain the whole configurational part of the 
probability distribution function in the model system. 
The information closed loop realized by the PI controller 
(26) then allows one to determine completely the inter- 
action potential. Conversely, if we are dealing with more 
complex systems, that contain further degrees of freedom 
beside the position of the center of mass of the atoms, a 
ME solution is always possible, which will correspond 
to an effective potential. If, however, the complete tar- 
get potential is sought, then it is necessary to match the 
relevant degrees of freedom of the systems with further 
involvement of information; for instance the experimen- 
tal three body correlation function and the inclusion of 
higher order terms in the definition of K would be neces- 
sary if a three body interaction is present. As a final re- 
mark, we point out that this inversion technique has been 
discussed assuming that the RDF of the target system is 
given. However, since experimental data are expressed in 
term of the structure factor, a preliminary transforma- 
tion to the real space RDF has to be performed in order 
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to extract the interaction potential of a real system. This 
procedure may be hampered by the limited range of the 
structure factor or by the unsatisfactory fc-resolution so, 
again, the use of the ME methods could reveal a useful 
tool to overcome those problems in optimal way. 
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